Pharmacy Deserts

The purpose of this notebook is to identify pharmacy deserts in Harris county.

Purpose: Accessibility to a pharmacy is an important facilitator of overall health, and is critical to ensure proper utilization of medications and adequate delivery of healthcase services. Community pharmacies offer a wide spectrum of services including dispensing medications, patient counselling, screening tests, immunization services, wellness programs, and education programs. Unlike many other healthcare settings, community pharmacies offer greater accessibility, with their extended hours of operation, availability of home delivery of medications, and no need to schedule an appointment for services.

However, not everyone has equal access to community pharmacies. Similar to “food deserts”, coined by the USDA, there exist “pharmacy deserts”, or geographic areas which lack access to a nearby pharmacy and where pharmacy services are scarece or difficult to obtain. Prominent public health researcher Dima M. Qato has defined pharmacy deserts in urban areas as regions where there is no pharmacy within one mile.

Methods: A dataset of licensed pharmacies in Texas was obtained from the Texas State Board of Pharmacy website. For the purposes of this analysis, pharmacies were limited to “Community Pharmacies” with active licenses. Other types of pharmacies, including those that only provide mail-order prescriptions, were excluded.

Census blocks are colored by closest distance from the center of the block to a pharmacy. Distances of greater than 1 mile have been used to denote pharmacy deserts in urban areas, however, for vulnerable populations (elderly, ill, families with young children, or people below the poverty line), even smaller distances can prove to be prohibitive.

Next Steps: Quantification and analysis of pharmacy deserts in Harris County and associated demographic information from desert regions.

For a large interactive version of the map, click here

#Load all required libraries
library(rgeos)
library(rgdal)
library(raster)
library(leaflet)
library(leaflet.extras)
library(tidyverse)
library(sp)
library(mapview)
library(png)

library(sf)
library(tidyverse)

library(maptools)
library(here)
library(stringr)
library(janitor)
library(readr)

library(plotly)
library(dplyr)
library(htmlwidgets)
library(rgeocodio)
library(geosphere)
library(PBSmapping)

#Obtain pharmacy dataframe, geolocate pharmacies, and obtain census information for each pharmacy

setwd("~/pharmacy_deserts")
p <- read.csv(file="phydsk.csv", header=T, stringsAsFactors=F) # 8,406 records

# filter the dataset to keep the pharmacies we want
unique(p$PHY.TYPE)

# look at subsets of data
unknown <- p[p$PHY.TYPE == "Unknown", ] # remove
mail.service <- p[p$PHY.TYPE == "Mail Service", ] # remove
other <- p[p$PHY.TYPE == "Other", ] # keep

# remove subsets of data
p <- p[p$CLOSED.DOOR == "N", ] # no closed door pharmacies (down to 7,741 records)
p <- p[!p$PHY.TYPE == "Unknown", ] # no "unknown" types of pharmacies (down to 7,710 records)
p <- p[!p$PHY.TYPE == "Mail Service", ] # no "mail service" types of pharmacies (down to 7,709 records)
p <- p[p$STATE == "TX", ] # only locations in Texas (down to 7,110 records)
p <- p[!p$LIC_STATUS == "Revocation", ] # remove revoked licenses (down to 7,054 records)

# drop unnecessary columns from dataset
colnames(p)
p[,c(1:2,5,10,12:13,18,20,22:27,29:35,45:51)] <- list(NULL) 

#Only Harris County
p <- p[p$COUNTY == "HARRIS",]
p <- p[p$LIC_STATUS == "Active",]
p <- p[p$CLASS == "Community Pharmacy",]

# geocode
library(rgeocodio)
geo <- gio_batch_geocode(paste(p$ADDRESS1,p$CITY,p$STATE,p$ZIP), fields = "census", api_key = gio_auth())

geo$lat <- geo$response_results
geo$lon <- geo$response_results

geo$county_fips <- geo$response_results
geo$census_tract <- geo$response_results
geo$block_group <- geo$response_results

geo$county_fips <- gsub("^.*county_fips","",geo$county_fips)
geo$county_fips <- gsub("\\,.*$","",geo$county_fips)
geo$county_fips <- gsub("\\=","",geo$county_fips)
geo$county_fips <- gsub("c\\(\\\"","",geo$county_fips)
geo$county_fips <- gsub("\\\"","",geo$county_fips)
geo$county_fips <- as.numeric(as.character(geo$county_fips))

geo$census_tract <- gsub("^.*tract_code","",geo$census_tract)
geo$census_tract <- gsub("\\,.*$","",geo$census_tract)
geo$census_tract <- gsub("\\=","",geo$census_tract)
geo$census_tract <- gsub("c\\(\\\"","",geo$census_tract)
geo$census_tract <- gsub("\\\"","",geo$census_tract)
geo$census_tract <- as.numeric(as.character(geo$census_tract))

geo$block_group <- gsub("^.*block_group","",geo$block_group)
geo$block_group <- gsub("\\,.*$","",geo$block_group)
geo$block_group <- gsub("\\=","",geo$block_group)
geo$block_group <- gsub("c\\(\\\"","",geo$block_group)
geo$block_group <- gsub("\\\"","",geo$block_group)
geo$block_group <- as.numeric(as.character(geo$block_group))

geo$GIDBG <- paste(geo$county_fips, geo$census_tract, geo$block_group, sep = "")

geo$lon <- gsub("^.*lng","",geo$lon)
geo$lon <- gsub("\\,.*$","",geo$lon)
geo$lon <- gsub("\\=","",geo$lon)
geo$lon <- gsub("c\\(","",geo$lon)
geo$lon <- as.numeric(as.character(geo$lon))

geo$lat <- gsub("^.*lat","",geo$lat)
geo$lat <- gsub("\\,.*$","",geo$lat)
geo$lat <- gsub("\\=","",geo$lat)
geo$lat <- gsub("c\\(","",geo$lat)
geo$lat <- as.numeric(as.character(geo$lat))

geo[,2:15, 18:20] <- list(NULL)
colnames(geo) <- c("geocoding_address","lat","lon", "GIDBG")
full <- cbind(p,geo)

write.csv(full,"pharmacies_geocoded_withblock.csv",row.names=F) # write the file

#Obtain shape file for Harris County

#Set up projections
wgs84<-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")
utm14n<-CRS("+proj=utm +zone=14 +ellps=GRS80 +datum=NAD83 +units=m +no_defs +towgs84=0,0,0")

#Census planning data ######
#Download from here: https://www.census.gov/research/data/planning_database/2018/
dfbg <- read_csv(here("pdb2018bgv4_us.csv")) %>% 
  dplyr::select(GIDBG:`Block_group`, 
                contains("ACS_12_16"), -contains("pct"))

#Prepare census data
harris_bg <- dfbg %>% 
  filter(State=="48" & County=="201") %>% 
  mutate(
    GEOID = GIDBG
  )

#If you don't want to re-run the earlier code to geocode the pharmacies, just uncomment this block and read the file back
#full <- read.csv(file="pharmacies_geocoded_withblock.csv", header=T) # read it back in
#full$GIDBG <- as.character(full$GIDBG)

#Combine pharmacy and census data
pharms.df <- left_join(full, harris_bg, by = "GIDBG")

#Import census blocks
library(tigris)
options(tigris_use_cache = FALSE)
blockg_sp <- block_groups(state = '48', county = '201',refresh = TRUE)

blockg_sp <- spTransform(blockg_sp, utm14n)

blockg_sf <- st_as_sf(blockg_sp) %>% 
  left_join(., harris_bg, by = "GEOID") %>% 
  mutate(
    `% Renting` = Renter_Occp_HU_ACS_12_16 / Tot_Occp_Units_ACS_12_16,
    `% Ages 18-24` = Pop_18_24_ACS_12_16 / Tot_Population_ACS_12_16,
    `% HH female head, no husband` = Female_No_HB_ACS_12_16 / Rel_Family_HHD_ACS_12_16,
    `% Non-White` = 1 - (NH_White_alone_ACS_12_16 / Tot_Population_ACS_12_16),
    `% Under age 65` = 1 - (Pop_65plus_ACS_12_16 / Tot_Population_ACS_12_16),
    `% HH with child under 6` = Rel_Child_Under_6_ACS_12_16 / Rel_Family_HHD_ACS_12_16,
    `% Vacant units` = Tot_Vacant_Units_ACS_12_16 / Tot_Housing_Units_ACS_12_16,
    `% Lacking college degree` = 1 - (College_ACS_12_16 / Pop_25yrs_Over_ACS_12_16),
    `% Hispanic` = Hispanic_ACS_12_16 / Tot_Population_ACS_12_16,
    `% Single housing units` = Single_Unit_ACS_12_16 / Tot_Occp_Units_ACS_12_16,
    `% Below poverty` = Prs_Blw_Pov_Lev_ACS_12_16 / Pov_Univ_ACS_12_16, 
    `% in diff. housing unit 1yr ago` = Diff_HU_1yr_Ago_ACS_12_16 / Tot_Population_ACS_12_16,
    `% Ages 5-17` = Pop_5_17_ACS_12_16 / Tot_Population_ACS_12_16,
    `% Black` = NH_Blk_alone_ACS_12_16 / Tot_Population_ACS_12_16,
    `% HH single person` = Sngl_Prns_HHD_ACS_12_16 / Rel_Family_HHD_ACS_12_16,
    `% Lacking HS degree` = Not_HS_Grad_ACS_12_16 / Pop_25yrs_Over_ACS_12_16,
    `% Crowded housing units` = Crowd_Occp_U_ACS_12_16 / Tot_Occp_Units_ACS_12_16,
    `% Without phone service` = Occp_U_NO_PH_SRVC_ACS_12_16 / Tot_Occp_Units_ACS_12_16,
    `% Over age 65` = Pop_65plus_ACS_12_16 / Tot_Population_ACS_12_16,
    
    high_elderly = ifelse(`% Over age 65`>=quantile(`% Over age 65`, probs=.75, na.rm=T) , 1, 0),
    high_below_poverty = ifelse(`% Below poverty`>=quantile(`% Below poverty`, probs=.75, na.rm=T) , 1, 0),
    high_rentals = ifelse(`% Renting`>=quantile(`% Renting`, probs=.75, na.rm=T) , 1, 0),
    high_kids6 = ifelse(`% HH with child under 6`>=quantile(`% HH with child under 6`, probs=.75, na.rm=T) , 1, 0),
    high_hispanic = ifelse(`% Hispanic`>=quantile(`% Hispanic`, probs=.75, na.rm=T) , 1, 0),
    high_black = ifelse(`% Black`>=quantile(`% Black`, probs=.75, na.rm=T) , 1, 0)
    
  )

#Convert data back to SF and proper projection for leaflet
mapdf <- spTransform(as(blockg_sf, "Spatial"), wgs84)

#Calculate centroids of each block group
harris.ps <- SpatialPolygons2PolySet(mapdf)
harris.centroids <- calcCentroid(harris.ps, rollup=1)
mapdf$CentroidX <- harris.centroids$X
mapdf$CentroidY <- harris.centroids$Y

#Determining the minimum distance from each block group centroid to pharmacy
#Note: This takes a very long time to run
for (i in 1:length(mapdf$GIDBG)){
  cents <- c(mapdf$CentroidX[i], mapdf$CentroidY[i])
  temp <- data.frame(Distance = integer())
  
  for (j in 1:length(pharms.df$lat)){
    temp <- rbind(temp, distVincentyEllipsoid(c(pharms.df$lon[j],pharms.df$lat[j]), cents))
  }
  mapdf$MinDistToPharm_Miles[i] <- min(temp)*0.00062137
}

#Create bins of minimum pharmacy distances for mapping
mapdf$pharm.bins <- ifelse(mapdf$MinDistToPharm_Miles <= 0.25, "0.00-0.25mi",
                           ifelse(mapdf$MinDistToPharm_Miles <= 0.5, "0.26-0.50mi",
                                  ifelse(mapdf$MinDistToPharm_Miles <= 1, "0.51-1.00mi",
                                         ifelse(mapdf$MinDistToPharm_Miles <= 2, "1.01-2.00mi",
                                                ifelse(mapdf$MinDistToPharm_Miles >2, "2.00mi+",NA)))))

rental <- subset(mapdf, mapdf@data$high_rentals==1)
kids <- subset(mapdf, mapdf@data$high_kids6==1)
hispanic <- subset(mapdf, mapdf@data$high_hispanic==1)
black <- subset(mapdf, mapdf@data$high_black==1)
elderly <- subset(mapdf, mapdf@data$high_elderly==1) 
poverty <- subset(mapdf, mapdf@data$high_below_poverty==1)

#Make map

pal <- colorFactor(c("#094d8c", "#2A92F2", "#ffaa5a","#FF7F05", "#af5500", "gray" ), 
                   domain = c("NA", "0.00-0.25mi", "0.26-0.50mi", "0.51-1.00mi", "1.01-2.00mi", "2.00mi+"))


#Popup box
popup<-paste0("<b>Block Group</b>: #", mapdf@data$GEOID, "<br/>",
              #"<b>% Renters</b>: ",  round(mapdf@data$X..Renting*100), "%<br/>",
              "<b>% Over 65 yrs</b>: ", round(mapdf@data$X..Over.age.65*100), "%<br/>",
              "<b>% HH with children <6</b>: ", round(mapdf@data$X..HH.with.child.under.6*100), "%<br/>",
              "<b>% Hispanic</b>: ", round(mapdf@data$X..Hispanic*100), "%<br/>",
              "<b>% Black</b>: ", round(mapdf@data$X..Black*100), "%<br/>",
              "<b>% Below Poverty</b>: ", round(mapdf@data$X..Below.poverty*100), "%"
)


map <- leaflet() %>%
  addProviderTiles(providers$OpenStreetMap) %>% 
  addPolygons(
    group = "All block groups",
    data = mapdf,
    fillColor = ~pal(mapdf@data$pharm.bins),
    popup = popup,
    weight = 0.5,
    opacity = 1,
    color = 'white',
    fillOpacity = 0.8
  ) %>% 
  addPolygons(
    group = "High % Over 65yrs",
    data = elderly,
    fillColor = ~pal(elderly@data$pharm.bins),
    popup = popup,
    weight = 0.5,
    opacity = 1,
    color = 'white',
    stroke = 0.1,
    fillOpacity = 0.8
  ) %>% 
  addPolygons(
    group = "High % HH kids <6",
    data = kids,
    fillColor = ~pal(kids@data$pharm.bins),
    popup = popup,
    weight = 0.5,
    opacity = 1,
    color = 'white',
    stroke = 0.1,
    fillOpacity = 0.8
  ) %>% 
  addPolygons(
    group = "High % Hispanic",
    data = hispanic,
    fillColor = ~pal(hispanic@data$pharm.bins),
    popup = popup,
    weight = 0.5,
    opacity = 1,
    color = 'white',
    stroke = 0.1,
    fillOpacity = 0.8
  ) %>% 
  addPolygons(
    group = "High % Black",
    data = black,
    fillColor = ~pal(black@data$pharm.bins),
    popup = popup,
    weight = 0.5,
    opacity = 1,
    color = 'white',
    stroke = 0.1,
    fillOpacity = 0.8
  ) %>%
    addPolygons(
    group = "High % Below Poverty",
    data = poverty,
    fillColor = ~pal(poverty@data$pharm.bins),
    popup = popup,
    weight = 0.5,
    opacity = 1,
    color = 'white',
    stroke = 0.1,
    fillOpacity = 0.8
  ) %>% 
  addMarkers(group = "Pharmacies", data=pharms.df, lat=~lat, lng=~lon, popup = ~as.character(pharms.df$PHARMACY_NAME)) %>%
  addLegend("bottomright", pal = pal,
            values = mapdf@data$pharm.bins,
            title = "Closest Pharmacy Distance<br/>by Block Group",
            opacity = 0.8) %>%
  addLogo(img = "http://januaryadvisors.com/wp-content/uploads/2013/05/ja-og.png", 
          position = c("bottomleft"), offset.x = 10, offset.y = 10, width = 90, height = 90, alpha = 0.9) %>% 
  addLayersControl(
    baseGroups = c("All block groups", "High % Over 65yrs", "High % HH kids <6", 
                   "High % Hispanic", "High % Black", "High % Below Poverty"),
    overlayGroups = c("Pharmacies"),
    options = layersControlOptions(collapsed = FALSE)
  )

# Export Map #####
htmlwidgets::saveWidget(map, here("PharmacyDeserts.html"))

#Generate map
map